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NONPARAMETRIC  ESTIMATION  OF  THE  PROBABILITY  OF 
A  LONG  DELAY  IN  THE  M/G/1  QUEUE 

BY 

D.  P.  GAVER 
P.  A.  JACOBS 

OPERATIONS  RESEARCH  DEPARTMENT 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  939^3-5000 

1 .  Introduction 

The  application  of  probability  theory  to  a  wide  variety  of  congestion 
problems  has  been  well  catalogued  in  many  papers.  Most  of  the  elegant 
solutions  obtained  are  in  somewhat  implicit  form,  being  presented  as 
functional  equations,  or,  frequently,  as  integral  (Laplace)  transforms, 
generating  functions,  and  sometimes  as  combinations  of  the  above.  The 
results  obtained  naturally  appear  in  terms  of  component  distribution 
functions  and  stochastic  processes  (renewal,  Poisson  etc).  Only  rarely  are 
issues  addressed  that  arise  when  actual  data  is  to  be  used  as  a  basis  for 
inference  from  the  models;  however,  see  Cox  [1965]. 

In  this  paper  we  consider  the  nonparametric  estimation  of  the 
probability  of  a  long  customer  delay  in  an  M/G/1  system,  given  a  known 
Poisson  arrival  rate  X  and  observations  of  independent  service  times  from 
the  service  distribution,  presumed  unknown.  Although  the  approach  and 
results  are  given  concretely  for  the  M/G/1  system,  they  apply  more  widely. 

To  be  specific,  consider  a  single  server  system  approached  by  a 
stationary  Poisson  (A)  traffic  with  A  known.  Service  times,  X  are 
independent  identically  distributed  with  unknown  distribution  function  F^ ; 
assume  AE[X]  -  p  <  1.  Let  observations  of  the  service  times  be  all  that  is 


known  about  F;  denote  these  by  x1  ,  x^,  . ..,  x^.  The  objective  is  to  supply 
nonparametric  estimates,  and  error  assessments  thereof,  of  the  probability 
of  a  long  delay  experienced  by  an  arriving  customer. 

It  is  well  known  that  if  W(t)  is  the  virtual  waiting  time  in  the  M/G/1 
queue  and  p  <  1 ,  then  the  moment  generating  function 


where  A(s)  is  the  moment  generating  function  of  a  distribution  H.  If  A(s) 
exists  for  s  <  sQ,  sQ  >  0,  then  there  will  be  a  smallest  real  zero  s  »  <  >  0 
of  the  denominator  of  (1.1)  which  can  be  used  to  show  that 

P{W  >  w}  -  D(k)  e-KW,  w  -*•  ®.  (1.2) 

We  will  always  assume  this  is  the  case. 

One  way  of  establishing  (1.2)  is  to  introduce 

3W  00 

¥(s)  -  — — —  -  |  P{W  >  w}  eSWdw  (1.3) 

0 

into  (1.1)  and  to  rewrite  in  the  form 

f(s)  -  pIMlLlJJ  +  pA(3)T(a)  (1.4) 
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which  is  equivalent  to  the  terminating  renewal  equation,  see  Feller  [1966] 
p.  362, 


w 

FyCw)  «  P{ W  >  w}  -  pH(w )  +  p  |  P{W  >  w-x)  H(dx)  (1.5) 

0 


where 


H(w) 


r  P  {  X  >  y} 

J  E[X] 

0 


dy. 


(1.6) 


—  ft  —  QW 

Following  Feller,  introduce  Fw  (w)  «  F^(w)  e  ’  ®  real  and  positive,  into 
(1.5)  to  obtain 


w 

Fw#(w)  -  p  H#(w)  ♦  J  F^(w-x)  pe0XH(dx). 

0 


Choose  9  =»  <  so  that 


CO  CO 

j  pe0XH(dx)  -  j  H#(dx)  -  1 
0  0 


(1.7) 


(1.8) 


yielding  a  standard  renewal  equation  for  F^. 
it  follows  that 


lim  Fy(w)  -  lim  F  ( w )e<W  - 

W  W  m(<) 

w-*®  w-*"=° 


From  the  key  renewal  theorem, 


(1.9) 


Availability  Codes 

Avail  and/or 

Qi  jt 

Special 

to 

X»/v 


too 


□  □ 


where 


c(<) 

and 


p  j  e**  H(x)dx 
0 


m(<)  -  p  J  xe<X  H(dx). 
0 

Summarizing 


P{W  >  t} 


c(<) 

m(<) 


~<t 


where  k  is  the  positive  solution  to  the  equation 


A0  1 C  4>(  0)  -  1]  -  1 


with 

4>(e)  -  E[e0X] . 


(1.10) 


(1.11) 


(1.12) 


(1.13) 


(1.14) 


In  section  2,  a  nonparametric  estimate,  <,  of  k  is  studied  which  solves 
equation  (1.13)  with  the  empirical  moment  generating  function  replacing  <f>. 
This  estimate  is  related  to  an  estimate  studied  by  Stigler  [1971]  in  the 
context  of  estimating  the  probability  of  extinction  for  branching  processes, 
ic  is  shown  to  be  a  consistent  estimate  of  k  having  a  distribution  in  the 
domain  of  attraction  of  a  stable  law.  Under  certain  conditions,  a  central 
limit  theorem  for  <  can  be  obtained. 
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In  section  3  a  non parametric  estimate  of 


iv 


y-' 

th 


-  c(<)  -tct 

p<t)  *  STkT  e 


(1.15) 


is  given. 

In  Section  4,  some  results  of  simulation  studies  of  the  estimates  of  < 
and  p(t)  are  presented. 

2.  Nonparametric  Estimation  of  the  Exponent  <  of  the  Probability  of  a 
Long  Wait 

Assume  for  the  remainder  of  the  paper  that  A»1  is  known.  Let 
x1 ,  ....  be  the  observations  of  the  service  times.  A  non-parametric 
estimate  of  the  moment  generating  function  of  the  service  time  distribution 
is 


.  n  ex 
<*>(  9)  =  —  I  e 
1  i»1 


(2.1) 


The  sample  equivalent  to  equation  (1.13)  is  thus 


i  -  e  1  [<f)(e)-i ]. 


(2.2) 


At  0  •  0,  the  RHS  of  (2.2)  is  x  which  is  less  than  1  if 


X1  +  X2  +  +  Xn 


<  1 ;  the  data  can  only  be  analyzed  for  a  stationary 


model  under  this  assumption,  which  will  be  made  in  what  follows.  Further, 
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L<V/\ 


(2.2)  has  a  unique  solution,  k,  which  is  an  estimate  of  <.  A  three-term 
Taylor's  expansion  of  the  RHS  of  (2.2)  about  8-0  gives  an  initial  estimate 


<H  -  2[ 1  -  x]  (x7)'1  (2.3) 

k  1  k  k  th 

where  x  -  —  [x,  ♦  ...  +  x  ],  the  k  sample  moment.  Since 
n  i  n 

^  1  A  ^  j  A  A 

0-1  L <t>( G )  -  1]  Z  x  +  -  8  x2,  <  is  an  upper  bound  for  k.  Equation  (2.2)  can 

<L  H 

be  solved  via  search  or  Newton-Raphson  iteration  to  obtain  the  estimate  c. 

We  will  now  present  asymptotic  results  for  the  distribution  of  tc  as  the 

sample  size  n  -*•  ®.  By  assumption,  if  X  is  a  random  variable  having  the 

B'X 

service  time  distribution,  then  E[e  ]  <  «  for  some  B'  >  k  >  0.  Thus, 
there  is  B  >  k  such  that  for  all  0  <  b  $  B 

E[XnebX]  <  (2.4) 

It  follows  from  the  monotonicity  of  the  RHS  of  (2.2)  that 

pu  >  si  -  .--t  -  , 

D  D  D 

-  PU(B)  -  4(B)  <  B[  1  -  .12]}.  (2.5) 

Since  $(0)  is  a  monotone  function  and  B  >  <,  [1  -  B  ^^(B)  -1]]  is  negative. 
Thus,  by  the  strong  law  of  large  numbers 
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(2.6) 


lim  P{<  >  B}  «  0. 
n-»® 


Let 


f  (6) 


0-1[J(e)  - 


1]. 


(2.7) 


Expand  f(<)  in  Taylor's  series  about  the  solution  <  of  (1.13).  Since 
f(<)  -  0, 

0  -  f(-c)  +  (<-<)  f'U)  ♦  1  (<-<)2  f"(0<)  (2.8) 


for  uome  6.  Thus 


-fUHf'U)  ♦ 


2  (< 


■<)  f '  '(0k)] 


(2.9) 


'  «-Y  ,,  y 

*  —  E)->fXe  +  e  -  1  ]  —  B  <  0 . 


(2.10) 


■ .  f  large  numbers 


it,  :  «• 


im  f ' (<)  -  8  <  0 


(2. 1 1  ) 


(2.12) 
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SI 


>*% 


m 


VS* 


with  probability  1.  It  follows  from  (2.6),  (2.9),  (2.11)  and  (2.12)  that  < 


converges  in  probability  to  k  as  n  -*■  ®.  Thus,  k  is  a  consistent  estimate  of 


If  E[e^<"<‘]  <  ®,  then 


Var[f(<)]  »  - - ^  VarLe*'*']  •  —  r2 

n  ic  n 


(2.13) 


By  the  central  limit  theorem 


—  /  n  { tc  -  <}  »  ^  /  n  [-f(<)]  [f  ’  ( <)  +  ^  (<~ic)  P'(0k)]  1  (2.19) 


is  asymptotically  standard  normal  as  n  ■*  *. 


2icX  <X 

If  E[e  ]  *  ”,  then  the  distribution  G  of  e  is  in  the  domain  of 


attraction  of  a  stable  law  with  index  1  <  a  <  2,  see  Feller  [1966].  Let 


{a^}  be  a  sequence  of  numbers  such  that 


a 

p  |  n  x!  G(dx)  1 
n  0 


as  n  -*■  °°. 


(2.15) 


The  normalized  random  variable 


^  .  n  k  X .  a  .  **  m 

—  [<  -  <3  -  —  [--  l  e  x][f’(<)  *  k<  -  <)  f "  ( e<)  ]~ 
3i  3  n  ,  d. 

n  n  i-1 


(2.16) 


.  .  .%  .*-  A  %  /•  .A  , 

“■**  •'»  i  J  *  ,  W  ,  ,  v  .  H  -  »  J  n  •»  _  >> 

»  »  »  »  -  V  -  It  _  »  »•  J %  y>»  _  »  „  *  -  »  T«  \ 


y  v 


is  in  the  domain  of  a  stable  law  with  index  a,  where  1  <  a  <  2  is  such  that 


y  ? 

j  x  G(dx) 
0 


2-cx 

y 


L(y)  where  L  is  a  slowly  varying  function. 


To  summarize,  we  have  the  following  result. 


PROPOSITION 


a.  k  is  a  consistent  estimate  of  <• 

2<X 

b.  If  E[e  ]  <  “,  then  the  distribution  of  <-<  is 

asymptotically  normal. 

2k:X 

c.  If  E[e  ]  =  ®,  then  <-<■  is  in  the  domain  of 
attraction  of  a  stable  law. 


We  will  suppose  for  the  remainder  of  this  section  that  the  service  time 


distribution  is  exponential  with  mean  ^  >  j.  In  this  case,  (since  by 


assumption) , 


<  -  u  *  1  ; 

Var[eKX]  »  u(u~1)2  ( 2- u )  '  ; 

and 

E[ f ' ( <) ]  =  -1. 

From  the  central  limit  theorem  (2.14) 
nVarU)  -  (^)2  =  u(2-u)  1  = 


(2.17) 

(2.18) 

(2.19) 

<*  1 
1  -IT  ' 


(2.20) 


Therefore,  if  g(x)  -  (1-x2)^/2  -  1  -  sin  V-x),  then 

nVar[g(<)]  -  [g' (<)]2nVar[<]  -  1. 

Thus,  g  is  an  asymptotically  variance-stabilizing  transform  of  k.  The 
simpler  related  transformation  ln(1+x)  is  used  in  the  simulations  of  < 
reported  in  Section  4. 

3.  A  Nonparametric  Estimate  of  P{W  >  t} 

The  asymptotic  analytic  result  is 

P(W  >  t}  -  e_<t  »  p(t)  (3.D 


as  t  -»  ® 


where 

00  CO 

c(<)  -  |  e<y  J  P{X  >  x}dx  dy 
0  y 

=  |  P{ Xedx}  ♦  ^  J  [e<X 
0  0 


(3-2) 


icx  -  P{  Xedx } 


and 


m(<)  -  j  xe<x  P ( X  >  x)  dx  (3-3) 

0 
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|J  tcxe<X  P{Xedx)  -  J  (e<X  -  1)  P(Xedx}} 

0  0 

00  00 

~  I  x2  P{Xedx}  +  ~  k  j  x3  PfXedx} 


(kx  -  1 )  R(x)  P(Xedx} 


where 


R(x)  »  [eKX  -  1  -  <x  - 


An  estimate  of  P{W  >  t}  is 


A  /  .  v  c ( <)  -«t 

p(t)  -  5-- —  e 


where  <  is  the  positive  solution  to  equation  (2.2); 


c(k)  -  ^  x2  +  5—  R1  (<) ; 


m(<)  »  ^  x2  ♦  ^  <  P  +  R^(<); 


-  -  n  <x .  -  ( <x .  ) 

R1(<)  -  -  I  Ce  1  -  1  -  <x  -  — 

i»1 


(3.4) 


(3.5) 


(3.6) 


(3.7) 


(3.8) 


W/VvVVvV>VvV 

-.VVaVaI 


m 


1 


(3.9) 


R.(k)  -  (-  R.(c))  +  ^  “  1  x. (e 
£  1-1  1 


n 


<x, 


KX, 


(<x.  ) 
i 


The  forms  of  c(0  and  m(<)  are  chosen  for  the  numerical  stability  of  the 

~  *  *  *  -i 

ratio  c(k)  m(<) 

In  the  remainder  of  this  section  we  will  assume  the  service  time 
distribution  is  exponential  with  mean  “  >  ^  and  x  <  1.  In  this  M/M/1  queue 
case,  it  is  well  known  that 


p(t)  -  exp{-(u~1)t}  -  (1  +  k)  1  e  Kt  (3.10) 

from  (2.16). 

We  will  now  motivate  a  transformation  of  p(t)  which  is  used  in  the 

Y 

simulation  studies.  Let  Y  -  ln(1+ic),  then  k  -  e  -  1  and 


p(t)  »  exp{-Y  -  [e^  -  1]t}. 


(3.11) 


from  (2.20).  Hence, 


nVar[h(Y)]  =  [-1  -e7t]2  [1  -  (e7  -  I)2]  1 

«  [-  1  -  t  -  (e7  -  1)t]2[1  -  (e7  -  I)2]"1.  (3.U) 

It  follows  from  the  definition  of  h(Y)  that 

nVar[h(Y) ]  *  [-(1+t)  +  h(Y)  +  Y]2{1  -  (e7  -  I)2}"1  .  (3.15) 

Since  0  <  k  -  e7  -  1  <  1  ,  0  <  Y  <  In  2.  Thus,  for  t  large 

nVar [h (Y) ]  -  [-h(Y)  +  t]2  (3-16) 

which  suggests  that  the  transformation  ln[t  -  In  p(t)J  will  tend  to 
stabilize  the  variance  of  p(t),  at  least  for  exponential  service  time 
distributions. 

4.  Simulation  Results 

In  this  section  results  are  presented  of  a  simulation  experiment  to 
study  small  sample  behavior  of  the  estimates  of  <  and  P{W  >  t}.  All 
simulations  were  carried  out  on  an  IBM  3033AP  computer  at  the  Naval 
Postgraduate  School  using  the  LLRANDOM  II  random  number  generating  package; 
Lewis  and  Uribe  [1981].  In  each  replication  50  exponential  service  times 
are  generated.  If  the  mean  of  the  service  times  is  larger  than  1,  equation 
(2.2)  has  no  positive  solution.  In  this  case,  another  50  service  times  are 
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generated.  The  estimates,  <  of  (2.2)  and  p(t)  of  equation  (3.5)  are 
computed  for  each  replication. 

Standard  errors  of  the  estimates  are  estimated  in  two  ways.  One  is  the 
appropriate  asymptotic  variance  expression.  The  other  is  the  jackknife 
procedure. 

The  jackknife  is  a  procedure  originally  introduced  by  Quenouille  for 
bias  reduction  [1956],  and  adapted  by  Tukey  [1958]  to  obtain  approximate 
confidence  intervals.  Suppose  interest  is  on  a  parameter  0  (e.g.  <  or  p(t)) 


that  is  estimated  by  0  using  a  complex  calculation  from  data  x^,...,x  .  The 


idea  is  that  of  assessing  variability  by  recomputing  0  after  removing 

independent  subgroups  of  data  of  equal  size,  and  then  using  the  recomputed  0 

values  to  estimate  a  variance  which  is  in  turn  applied  to  state  a  standard 

error  or  a  two-sided  confidence  interval  that  contains  the  true  0  with 

specified  confidence.  A  few  more  details  follow;  for  more,  see  Efron  [1982] 

and  his  more  recent  work,  or  Mosteller  and  Tukey  [1977].  The  actual 

calculation  involves  splitting  the  n  data  points  into  g  disjoint  groups  of 

size  m;  n-mg.  In  our  simulations  g  is  always  10.  Then  calculate 

j*1,2,...,g:  the  estimate  of  0  based  on  a  reduced  data  set  that  omits  the 
.  th 


j  group.  In  the  simulations,  the  first  group  is  the  first  five 
(unordered)  service  times;  the  second  group  is  the  second  five  service 
times,  etc. 

Now  Tukey  computes  pseudo-values 


yj  - 88  -  (8-°  6<-j) 
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which  are  treated  as  independent.  Tukey  recommends  referring  the  mean  of 


the  pseudo-values  y  to  Student's  t  with  g-1  degrees  of  freedom  to  obtain 
confidence  limits. 

In  the  jackknife  procedure  to  estimate  k,  it  is  sometimes  the  case  that 

a  positive  solution  of  equation  (2.2)  does  not  exist  for  the  data  set  that 

omits  the  subgroup  because  the  mean  of  this  reduced  data  set  is  larger 

than  1.  In  this  case,  k,  . .  is  set  equal  to  the  smallest  of  the  's  that 

(-J)  (1) 

could  be  computed  for  the  sample  considered.. 

Similiarly,  in  a  jackknife  procedure  to  estimate  p(t),  it  is  sometimes 

A  A 

the  case  that  either  <(_j)  does  not  exist  or  P(_j)(t)  exceeds  1  for  the 
reduced  data  set  that  omits  the  jth  subgroup.  In  this  case  P(_j)  i-3  set 
equal  to  the  largest  of  those  p^^'s  that  could  be  computed  and  were  less 
than  1 . 


4. 1  Results  of  a  Simulation  Experiment  to  Estimate  k 

In  this  subsection  results  are  given  of  a  simulation  experiment  to 
estimate  k.  For  each  replication,  three  different  estimates  of  <  were 

A 

computed.  Estimate  I,  is  computed  by  numerically  solving  (2.2)  by 
search  starting  with  the  initial  value  <u  of  (2.3).  Estimate  II,  <Tr 
is  obtained  by  jackknifing  k  using  ten  subgroups;  the  mean  of  the  pseudo¬ 
values  is  the  estimate.  Estimate  III  is  obtained  by  jackknifing  ln(l+<) 
using  ten  subgroups;  the  inverse  transform  of  the  mean  of  the  pseudo- values 

y 

e  -  1  is  used  as  the  estimate  of  k.  In  the  cases  of  Estimates  II  and  III, 


if  the  estimate  is  negative,  it  is  set  equal  to  0.  For  each  estimate,  the 
average  bias 


B 


1 

500 


500 

l 

i-1 


(*,  -tc) 


and  the  relative  mean  square  error 


Rel  MSE 


1 

500 


500  k. -k  2 

i  HH 

i-i 


are  computed  where 


IC  *  vi—  1 

in  the  case  of  exponential  service  times  with  mean  -  <  1. 

y 

Results  of  the  simulation  appear  in  Table  1.  All  of  the  estimates  have 
about  the  same  relative  mean  square  error.  Jackknife  estimates  II  and  III 

a 

have  smaller  bias  than  the  straightforward  estimate  I.  Jackknifing  k  itself 

rather  than  ln(x+1)  gives  the  smallest  bias.  As  -  increases  all  the 

y 

estimates  have  increased  relative  mean  square  error. 

A  simulation  study  was  conducted  to  compare  the  performance  of 
different  confidence  interval  procedures.  For  each  replication  three  80? 
confidence  intervals  were  constructed.  Interval  procedure  I  is  a  normal 
confidence  interval  which  uses  the  straightforward  estimate  of  <,  <  ,  as  the 
point  estimate;  the  estimate  of  the  variance  is  the  data  version  of  the 


asymptotic  variance  in  the  central  limit  theorem  (2.14); 


“2 

°CLT 


~  O  ~  _  O  1  “  <* 

K  B  — j-  l  (e  -  $(<))' 

i 


(H.1) 


where 


*(<)  -  -  I  e  (4.2) 

-  .  <X.  ICX. 

B  -  -  I  e  1  -  1  -  <x  e  1  (4.3) 

i 

with  <  «  <  In  this  case.  The  80%-point  of  the  normal  distribution  is  used 
to  construct  the  interval.  Limits  that  are  negative  are  set  equal  to  0. 

Confidence  interval  procedure  II  is  a  jackknife  confidence  interval 
which  jackknifes  <  and  uses  the  80%  point  of  the  student  t-distribution  with 
9  degrees  of  freedom.  Limits  that  are  negative  are  set  equal  to  0. 

Confidence  interval  procedure  III  i3  a  jackknife  confidence  interval 
wnich  jackknifes  ln(<+1)  and  uses  the  80%  point  of  the  student 
t-distribution  with  9  degrees  of  freedom  to  give  a  confidence  interval  for 
ln(<+1).  The  inverse  transformation  of  the  endpoints  of  the  interval  gives 
an  interval  for  <;  limits  that  are  negative  are  set  equal  to  0. 

Results  of  the  confidence  interval  simulation  appear  in  Table  2. 
Reported  are  the  number  of  500  intervals  that  cover  the  true  value  of 
<,  (C);  the  number  of  the  500  intervals  such  that  the  entire  interval  lies 
below  <,  (3);  and  the  number  of  the  500  intervals  such  that  the  entire 
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interval  lies  above  the  true  value,  (H).  The  average  length  of  the 
confidence  intervals  is  also  given. 

The  number  of  intervals  that  cover  the  true  value  of  <  for  procedure 

III  is  within  .80  ±  ( 1 .96)X  =7rr-  (  .2)  {  .8)  -  [.765,  .835].  All  but  one  case 

5UU 

of  the  normal  confidence  intervals  of  procedure  I  are  outside  this  range. 
All  but  one  case  using  confidence  interval  procedure  II  are  inside  this 
range.  The  average  width  of  confidence  intervals  for  procedures  II  and  III 
are  about  the  same.  Thus,  although  the  jackknife  estimate  tc  is  a  little 
more  biased  than  k  ,the  coverage  of  the  jackknife  confidence  interval  for 
appears  to  be  somewhat  better  than  for  jackknife  confidence  interval 


procedure  II. 


TABLE  1 

Bias  and  Mean  Relative  Square  Error 
for  Estimates  of  < 


Distribution 

Estimate 

Exponential 

1  -  0.6 
<  -  .6667 

Exponential 

1  »  o.7 
u 

<  *  .4286 

Exponential 

-  -  0.8 

u 

<  -  .2500 

Exponential 

-  -  0.9 

if  -  .1111 

B  Rel  MSE 

B  Rel  MSE 

B  Rel  MSE 

B  Rel  MSE 

(S.E) 

(S.E. ) 

(S.E) 

(S.E. ) 

(S.E) 

(S.E.) 

(S.E)  (S.E.) 

(-] 

(§) 

(-) 

(-) 

K 

K 

< 

I 

.1147 

.2350 

.0557 

.3292 

.0057 

.6026 

.0865 

2.400 

(.014) 

(.021 ) 

(.011) 

(.029) 

(.008) 

(.054) 

(.0067) 

(.2134) 

[.172] 

[.130] 

[.023] 

[.779] 

II 

.0219 

.2256 

.0067 

.3095 

.0102 

.5382 

.0469 

1  .785 

(.014) 

(.019) 

(.011) 

(.025) 

( .008) 

(.046) 

( .0063) 

(.161  ) 

[.033] 

[.016] 

[.041  ] 

[.422] 

III 

.0533 

.221  2 

.0160 

.3073 

.0268 

.5576 

.0608 

1.997 

(.014) 

(.019) 

(.011) 

(.027) 

( .008) 

( .049) 

( .008) 

(.177) 

[.080] 

[.037] 

[.1072] 

[.5472] 

Table  2 


Confidence  Intervals  for  < 


Procedure 

Distribution 

Coverage 

Width 

B 

% 

I 

(C) 

(%) 

H 

% 

B 

% 

II 

(C) 

(J) 

H 

% 

B 

% 

III 

(C) 

(%) 

H 

% 

I 

AV 

(S.E.) 

II 

AV 

(S.E. 

III 

AV 

(S.E. ) 

Exponential 

I  -  0.6 

45 

(363) 

92 

38 

(408) 

54 

29 

(395) 

76 

.6402 

.8089 

.7752 

ic  »  .6667 

9 

(72.6) 

18.4 

7.6 

(81 .6) 

10.8 

5.8 

(79) 

15.2 

( .005) 

( .01 24) 

(.0114) 

Exponential 

-  -  0.7 

\x 

49 

(380) 

71 

50 

(409) 

41 

42 

(399) 

59 

.5279 

.6255 

.6157 

<  »  . 4286 

9.8 

(76) 

14.2 

10 

(81.8) 

8.2 

8.4 

(79.8) 

11.8 

(.004) 

(.009) 

(.009) 

Exponential 

-  -  0.8 

11 

28 

(407) 

65 

38 

(420) 

42 

33 

(413) 

54 

.4441 

.4831 

.4903 

5.6 

(81 .4) 

13 

7.6 

(84) 

8.4 

6.6 

(82.6) 

10.8 

(.005) 

( .008) 

( .008) 

Exponential 

-  -  0.9 

JJ 

1 

(429) 

70 

54 

(408) 

38 

52 

(392) 

56 

.3733 

.3763 

•  3923 

<  -  .1111 

.002(85.8) 

14 

10.8 

(81.6) 

7.6 

10.4 

(78.4)11.2 

(.005) 

(.009) 

(.009) 

4 . 2  Simulation  Results  for  Estimating  p(t) 

In  this  subsection,  results  are  given  of  a  simulation  experiment  to 
estimate  p(t).  For  each  replication  three  estimates  of  p(t)  are  computed. 
Estimate  I,  p  (t)  is  computed  for  each  replication  using  formulas  (2.2), 
(3-5)  -  (3*9).  If  Pj(t)  exceeds  1,  it  is  set  equal  to  1. 
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There  are  (at  least  two)  possible  ways  to  implement  a  jackknife 

procedure  to  estimate  p(t).  Estimate  II  is  obtained  by  jackknifing  ln(«r+1); 

an  estimate  of  <  is  obtained  by  the  inverse  transformation  of  the  mean  of 

the  pseudo- values  y  , 

LK 


LK 

<n  -  e  -  1; 


Kjj  is  used  in  formulas  (3-5)  -  (3-9)  to  obtain  the  estimate  p  (t).  If 


p  (t)  exceeds  1,  it  is  set  equal  to  1. 


Estimate  III  is  obtained  by  jackknifing  ln[t  -  In  p ( t ) ] ;  if  p(t) 

.  th 


exceeds  1  f or  a  reduced  data  set  that  omits  the  j  subgroup,  the  estimate 
of  ln[t  -  In  p(t)]  for  that  reduced  data  set  is  put  equal  to  the  smallest 
estimate  that  could  be  computed  from  the  other  reduced  data  sets.  An 
estimate  of  p(t)  is  obtained  by  inverse  transformation  of  the  average  of  the 
pseudo-values  yLPR: 

-y, 


/ .  .  t  LPR 
PIII(t)  °  e  e 


If  p  ^(t)  >1,  it  is  replaced  by  1 


For  each  estimate,  the  average  bias 


B 


1 

500 


500  „ 

I  [p.  (t )  -  p(t )] 
i-1 


(«.«) 


and  the  relative  mean  square  error 


The  entries  for  estimate  II  in  Table  3  suggest  tr.at  ;  : 
then  computing  the  estimate  of  the  probability  using  tr.--  . 

of  <  increases  the  bias  and  relative  mean  square  --rr  :■  : 
p  (t)  the  strai ght-f orward  estimate.  The  entries  f  •  '  • 

the  jackknifing  In[t  -  In  p  ( t )  ]  gives  about  tne  ■- 
error  as  the  original  estimate  p  ( t )  but  incr-i 

A  simulation  experiment  was  conducted  tc  .  r '  ■ 
confidence  intervals  obtained  by  jackknifing 
simulation  replication,  the  average  and  van  an  -  :  ■  • 

computed; 


-  i_  y° 

yj  '  10  ^  y(j 


i 


10 


■  9  L,  (,(J)’y)' 


An  dot  confidence  interval  for  ln[t  -  In  p(t)]  is 


[LPL.LPV]  =  y  ±  (  1.383)  — - —  /"o* 

J  /  10  J 

where  1.383  is  the  80t  point  for  a  Student's  t-distribution  with  9  degrees 
of  freedom.  The  limits  of  the  interval  are  inversely  transformed  to  give  a 


confidence  interval  for  p(t). 


1^:1 


$ 


PL  »  exp{t  -  e 


P V  »  exp  (t  -  eLPIj} , 


If  a  confidence  limit  exceeds  1,  the  limit  is  set  equal  to  1. 

Results  of  the  confidence  interval  simulation  appear  in  Table  4. 
Reported  are  the  number  of  the  500  intervals  that  cover  the  true  value 
p(t)  (C);  the  number  of  the  500  intervals  for  which  the  upper  limit  PU  is 
below  p(t)  (B);  and  the  number  of  the  500  intervals  such  that  the  lower 
limit  PL  is  above  the  true  value  (H).  The  average  width  of  the  confidence 
interval  is  also  given.. 


Table  4. 

80$  Confidence  Intervals  for  p(t) 


Coverage 


Exponential 

i  T-3 

p ( 3 )  -  .3779 


B 

C 

H 

(l) 

(?) 

(?) 

'  56 

41  5 

29 

(11.2) 

(83) 

(5.8) 

52 

408 

40 

(10. U) 

(81 .6) 

(8) 

54 

412 

34 

(10.8) 

(82.4) 

(6.8) 

Average  Width 


(standard  error) 


Exponential 

^  -.6,  T-1  56  415  29  .325 

p ( 1 )  »  .3081  (11.2)  (83)  (5.8)  (.006) 


Exponential 

-  -.7,  T»2 

u 


The  coverage  of  the  confidence  intervals  is  within 
.80  ±  (1.96)  /(  -  .80  ±  .035  »  [.765,  .835].  The  width  of  the 


interval  increases  as  -  increases. 

y 


J.3  Simulation  Results  for  Estimating  p(t)  for  Mixed  Exponential  Service 


Times 

In  this  subsection,  results  are  given  of  a  simulation  experiment  to 
estimate  p(t)  -  P { W> t }  in  the  case  of  mixed  exponential  service  times. 


each  replication  50  random  numbers  are  generated  from  a  mixed  exponential 
distribution  with 


P(S>t) 


+ 


t>0. 


The  estimate  ln[t  -  lnp(t)]  is  jackknifed  with  ten  subgroups  as  before  where 
p(t)  is  given  by  (3.5)  and  80$  confidence  intervals  are  constructed.  Each 
simulation  has  500  replications. 

The  asymptotic  distribution  of  the  virtual  waiting  time,  W,  can  be 
found  by  inverting  the  moment  generating  function  (1.1);  it  is  a  mixture  of 
two  exponential  random  variables. 

Table  5  reports  results  of  three  simulations. 


Case  I : 


1 

E(S) 


.30 


“  0.6, 


P{ W>T}  - 


T 


-  3404  ; 


1  , 


Case  1 1 : 


ui 

E(S) 


•  35 


0.7, 


-  -  1.05,  T 
U2 

P { W> T)  -  .3^59; 


2, 


Case  III : 


-  .40 


-  -  1.20,  T  -  3, 

y2 


E(S)  -  0.8,  P{ W>T} 


.^333  ; 


For  each  simulation,  the  mean  bias  and  relative  mean  square  error  (4.4) 
and  (4.5)  are  computed  with  p(t)  =  P { W> t } .  Confidence  interval  coverage  and 
average  widths  are  also  computed. 

Both  cases  II  and  III  each  had  one  replication  for  which  using  all  the 
data  to  compute  the  estimate  p(t)  of  (3.5)  resulted  in  a  value  larger  than 
1;  these  two  replications  are  not  counted  in  the  summary  statistics  in  Table 

5. 

Comparison  of  the  mean  biases  and  mean  relative  square  errors  of  Table 
5  with  those  of  estimate  III  in  Table  3  shows  that  they  are  about  the  same 
for  both  service  time  distributions.  The  coverage  of  the  confidence 
intervals  in  Table  5  is  again  within  [.765,  .835].  However,  the  average 
length  of  the  confidence  intervals  for  the  mixed  exponential  cases  are 
larger  than  those  reported  in  Table  4  for  the  exponential  distributions. 


Table  5 


Results  for  Estimates  of  P { W> t } 
for  an  M/C/1  Queue  with  Mixed 
Exponential  Service  Times 


Distribution 

Case 

Estimate 

Bias  Rel  MSE 

(S.E)  (S.E.) 

Conf idence 
Coverage 

B  (C)  H 

}  (})  } 

Intervals 

Width 

Average 

(S.E.) 

I 

.024 

.172 

58 

(407) 

35 

.392 

( .006) 

(.01 3) 

11.6 

(81.4) 

7 

( .007) 

II 

.040 

.384 

63 

(408) 

28 

.500 

(.009) 

( .028) 

12.6 

(81  .4) 

5.6 

(  .010) 

III 

.022 

.350 

54 

(411) 

34 

.637 

(.011) 

(.021 ) 

10.8 

(82.4) 

6.8 

(.011) 

4.4  Simulation  Results  for  Estimating  P { W> t }  for  Gamma  Distributed  Service 
Times 

In  this  subsection,  results  are  given  of  a  simulation  experiment  to 

estimate  p(t)  =*  P{W>t}  in  the  case  of  gamma  service  times.  Each  experiment 

has  500  replications.  For  each  replication,  50  service  times  are  generated 

having  the  distribution  of  the  sum  of  two  exponential  random  variables  each 
1 

having  mean  -.  The  estimate  ln[t  -  In  p(t)]  is  jackknifed  with  ten 
subgroups  as  before  where  p(t)  is  given  by  (3.5)  and  80}  confidence 


intervals  are  constructed. 


The  asymptotic  distribution  or'  the  virtual  waiting  time,  W,  can  be 
found  by  inverting  the  moment  generating  function  (1.1);  it  is  a  mixture  of 
two  exponential  random  variables. 

Table  6  reports  results  of  three  simulations. 

Case  I :  -  -  . 30  T  «  1 , 

u 

ECS)  -  0.6,  P{W>T}  -  .2508; 

Case  II:  --.35  T  -  2, 

u 

ECS)  -  0.7,  P{ W>T}  -  .2232; 

Case  III :  -  -  .U0  T  »  3, 

u 

EC S)  -  0.8,  P { W> T }  -  .2950; 


VV 


Table  6 


Results  for  Estimates  of  P(W>t} 
for  an  M/G/1  Queue  with  Gamma 
Service  Times 


Distri bution 

“i 

Estimate 

Confidence 

Coverage 

Intervals 

Width 

Case 

Bias 

Rel  MSE 

B 

(C) 

H 

Average 

(S.E) 

(S.E. ) 

1 

(*) 

% 

(S.E.) 

I 

.007 

.110 

61 

(901  ) 

38 

.227 

( .0014) 

(  .007) 

1  2.2 

(80.2) 

7.6 

( .003) 

II 

.0167 

.299 

61 

(397) 

92 

.  328 

( .005) 

( .026) 

12.2 

(79.9) 

8.9 

(  .007) 

III 

.0251 

.390 

62 

(  907) 

31 

.990 

(.008) 

( .039) 

12.9 

(81.9) 

6.2 

( .010) 

The  mean  bias  and  relative  mean  square  error  are  about  the  same  as  for 
the  exponential  and  mixed  exponential  service  time  distributions.  The 
confidence  interval  coverage  is  once  again  within  [.765.  .835].  The  average 
width  of  the  confidence  intervals  is  smaller  than  the  widths  in  the 
exponential  and  mixed  exponential  cases.  This  is  to  be  expected  since  the 
gamma  distribution  has  a  shorter  tail  than  the  exponential  or  mixed 
exponential  distribution. 
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